In [103]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Immersed hot rectangular body:
In this part, we intend to place the body (rectangular or square) in our domain. What changes do you think need to be considered? Consider the hot rectangle that is located in the middle of the cavity however the vertical position of its center can be varied from 0.25 of the cavity height to 0.75 . The length of the immersed body is twice of its height. We run our code for three different Rayleigh number while the rectangle is in different positions.
Changes: We do not need to create nodes for the rectangle instead we use the lattice nodes which have already been created. We only add some thermal and fluid boundary conditions.
At first the location of left, right, top and bottom edges are determined. Find it in the code. Then we apply a constant temperature for all the sides in the section where thermal distribution function is solving. (Do you remember how we apply hot temperature for the left wall in our previous code?) . You can use it here for all the walls. All of the sides should be considered as solid and no-slip BCs need to be applied. (Do you remember how we apply no-slip for the walls in natural convection or liddriven cases?) You can exactly use it for all the sides where fluid distribution function (f) is solving.
Further experience:
Try to change the obstacle’s length. Change its position in the cavity. Try to change the hot obstacle to cold one. What happened?
Question:
If you notice we determine the number of lattice nodes in a way that the obstacle’s edges are located right on the lattice nodes which have been already created. Now the question is that how can we create a circle in the domain? One way is that you just make a condition for this circle in the code and bounce back conditions only be considered for the nodes which are in the circle. But it is a rough estimation and undoubtedly leads to errors and on the other side we need to increase the lattice nodes to decrease this error. In that case the lattice nodes are not directly located on the circle’s surface. Please look at the related papers to find more about this issue.
In [104]:
nx=48 # the number of nodes in x direction lattice direction
ny=48 # the number of nodes in y direction lattice direction
pr=0.71
visco=0.02 #*ny
dt=1.0 # temperature difference between left and right wall
alpha=visco/pr # heat diffusion coefficient
Tw=1.0 # left wall temperature
rhoo=1.0
D=3.0 # the dimension of the problem
omega=1./(0.5+(visco*D)) #Chapman-Enskog dt=1 and dx=1 relaxation flow
omegat=1./(0.5+(alpha*D)) # relaxation temperature
mstep=50000 # the number of time step
k=9 # k=0,1,2,3,4,5,6,8,9
In [105]:
w=numpy.ones(k) # witghting factor
cx=numpy.ones(k)
cy=numpy.ones(k)
rho=numpy.ones((ny+1,nx+1) ) # density matrix
u=numpy.ones((ny+1,nx+1) )
v=numpy.ones((ny+1,nx+1) )
f= numpy.ones((k, ny+1,nx+1)) # distribution function
feq= numpy.ones((k, ny+1,nx+1))
g= numpy.ones((k, ny+1,nx+1)) # distribution function for temperature
geq= numpy.ones((k, ny+1,nx+1))
T=numpy.ones((ny+1,nx+1) )
t1=numpy.ones((ny+1,nx+1) )
t2=numpy.ones((ny+1,nx+1) )
force=numpy.ones((ny+1,nx+1) )
In [106]:
strf=numpy.ones((ny+1,nx+1) ) # streamfunction
In [107]:
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0
cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0
In [108]:
##================== Initial value
rho[0:ny+1,0:nx+1]=rhoo
v[0:ny+1,0:nx+1]=0.0
u[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=Tw
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*rho[j,i]
g[l,j,i]=w[l]*T[j,i]
Tref=0.5
In [109]:
## Main loop : comprised two parts :collision and streaming
##=====================
def function (Ra,yc) :
#************ Obstacle lengths *************
xc=nx/2
#yc=ny/2 # ny/4,ny/2,3ny/4
yc=yc
dx=nx/12 # fixed = 1/10 of the cavity length
dy=ny/24 # fixed=1/20 of the cavity length
x1=xc-dx
x2=xc+dx
y1=yc-dy
y2=yc+dy
gb=(Ra*visco*alpha)/(ny**3)
for n in range(mstep) :
if (n%10000==0):
print n
# collision process fluid flow
t1[0:ny+1,0:nx+1]=u[0:ny+1,0:nx+1]*u[0:ny+1,0:nx+1]+v[0:ny+1,0:nx+1]*v[0:ny+1,0:nx+1]
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
force[0:ny+1,0:nx+1]=(3.*w[l]*gb*(T[0:ny+1,0:nx+1]-Tref)*cy[l]*rho[0:ny+1,0:nx+1]*cos(alpha))+\
(3.*w[l]*gb*(T[0:ny+1,0:nx+1]-Tref)*cy[l]*rho[0:ny+1,0:nx+1]*sin(alpha))
feq[l,0:ny+1,0:nx+1]=w[l]*rho[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1]+4.5*t2[0:ny+1,0:nx+1]**2-1.5*t1[0:ny+1,0:nx+1] )
force[0,0:nx+1]=0.
force[ny,0:nx+1]=0.
force[0:ny+1,0]=0.
force[0:ny+1,nx]=0.
f[l,0:ny+1,0:nx+1]=(1.-omega)*f[l,0:ny+1,0:nx+1]+omega*feq[l,0:ny+1,0:nx+1]+force[0:ny+1,0:nx+1]
f[1,0:ny+1,nx:0:-1]=f[1,0:ny+1,nx-1::-1]
f[3,0:ny+1,0:nx]=f[3,0:ny+1,1:nx+1]
f[2,ny:0:-1,0:nx+1]=f[2,ny-1::-1,0:nx+1]
f[5,ny:0:-1,nx:0:-1]=f[5,ny-1::-1,nx-1::-1]
f[6,ny:0:-1,0:nx]=f[6,ny-1::-1,1:nx+1]
f[4,0:ny,0:nx+1]=f[4,1:ny+1,0:nx+1]
f[8,0:ny,nx:0:-1]=f[8,1:ny+1,nx-1::-1]
f[7,0:ny,0:nx]=f[7,1:ny+1,1:nx+1]
## =============================
#left Boundary- Bounce back
f[1,0:ny+1,0]=f[3,0:ny+1,0]
f[5,0:ny+1,0]=f[7,0:ny+1,0]
f[8,0:ny+1,0]=f[6,0:ny+1,0]
#right Boundary-bounce back
f[3,0:ny+1,nx]=f[1,0:ny+1,nx]
f[7,0:ny+1,nx]=f[5,0:ny+1,nx]
f[6,0:ny+1,nx]=f[8,0:ny+1,nx]
# top Boundary -bounce back
f[4,ny,0:nx+1]=f[2,ny,0:nx+1]
f[7,ny,0:nx+1]=f[5,ny,0:nx+1]
f[8,ny,0:nx+1]=f[6,ny,0:nx+1]
# bottom Boundary -bunce back
f[2,0,0:nx+1]=f[4,0,0:nx+1]
f[5,0,0:nx+1]=f[7,0,0:nx+1]
f[6,0,0:nx+1]=f[8,0,0:nx+1]
#============= obstacle boundary =====
#right Boundary-bounce back
f[3,y1:y2+1,x2]=f[1,y1:y2+1,x2]
f[7,y1:y2+1,x2]=f[5,y1:y2+1,x2]
f[6,y1:y2+1,x2]=f[8,y1:y2+1,x2]
#left Boundary- Bounce back
f[1,y1:y2+1,x1]=f[3,y1:y2+1,x1]
f[5,y1:y2+1,x1]=f[7,y1:y2+1,x1]
f[8,y1:y2+1,x1]=f[6,y1:y2+1,x1]
# top Boundary -bounce back
f[4,y2,x1:x2+1]=f[2,y2,x1:x2+1]
f[7,y2,x1:x2+1]=f[5,y2,x1:x2+1]
f[8,y2,x1:x2+1]=f[6,y2,x1:x2+1]
# bottom Boundary -bunce back
f[2,y1,x1:x2+1]=f[4,y1,x1:x2+1]
f[5,y1,x1:x2+1]=f[7,y1,x1:x2+1]
f[6,y1,x1:x2+1]=f[8,y1,x1:x2+1]
for i in range(nx+1):
for j in range(ny+1):
sumr=0.0
for l in range (k):
sumr=sumr+f[l,j,i]
rho[j,i]=sumr
for i in range(1,nx):
for j in range(1,ny):
usum=0.0
vsum=0.0
for l in range (k):
usum=usum+f[l,j,i]*cx[l]
vsum=vsum+f[l,j,i]*cy[l]
u[j,i]=usum/rho[j,i]
v[j,i]=vsum/rho[j,i]
u[y1:y2+1,x1:x2+1]=0.0
v[y1:y2+1,x1:x2+1]=0.0
u[:,0]=0.
u[:,nx]=0.
u[0,:]=0.
u[ny,:]=0.
v[:,0]=0.
v[:,nx]=0.
v[0,:]=0.
v[ny,:]=0.
# collision process Temperature
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
geq[l,0:ny+1,0:nx+1]=w[l]*T[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1] )
g[l,0:ny+1,0:nx+1]=(1.-omegat)*g[l,0:ny+1,0:nx+1]+omegat*geq[l,0:ny+1,0:nx+1]
#streaming process
# ==========================
# streaming Temperature
g[1,0:ny+1,nx:0:-1]=g[1,0:ny+1,nx-1::-1]
g[3,0:ny+1,0:nx]=g[3,0:ny+1,1:nx+1]
g[2,ny:0:-1,0:nx+1]=g[2,ny-1::-1,0:nx+1]
g[5,ny:0:-1,nx:0:-1]=g[5,ny-1::-1,nx-1::-1]
g[6,ny:0:-1,0:nx]=g[6,ny-1::-1,1:nx+1]
g[4,0:ny,0:nx+1]=g[4,1:ny+1,0:nx+1]
g[8,0:ny,nx:0:-1]=g[8,1:ny+1,nx-1::-1]
g[7,0:ny,0:nx]=g[7,1:ny+1,1:nx+1]
## Boundary condition
##============================ Temp boundary
#
#left Boundary
g[1,1:ny,0]=-g[3,1:ny,0]
g[5,1:ny,0]=-g[7,1:ny,0]
g[8,1:ny,0]=-g[6,1:ny,0]
#right Boundary
g[3,1:ny,nx]=-g[1,1:ny,nx]
g[7,1:ny,nx]=-g[5,1:ny,nx]
g[6,1:ny,nx]=-g[8,1:ny,nx]
# top Boundary
for l in range(k):
g[l,ny,0:nx+1]=g[l,ny-1,0:nx+1]
# bottom Boundary
for l in range(k):
g[l,0,0:nx+1]=g[l,1,0:nx+1]
#
#
##========= obstacle temperature boundary ==============
#left Boundary
g[1,y1:y2+1,x1]=( Tw*(w[1]+w[3]) )-g[3,y1:y2+1,x1]
g[5,y1:y2+1,x1]=( Tw*(w[5]+w[7]) )-g[7,y1:y2+1,x1]
g[8,y1:y2+1,x1]=( Tw*(w[8]+w[6]) )-g[6,y1:y2+1,x1]
#right Boundary
g[3,y1:y2+1,x2]=( Tw*(w[1]+w[3]) )-g[1,y1:y2+1,x2]
g[7,y1:y2+1,x2]=( Tw*(w[5]+w[7]) )-g[5,y1:y2+1,x2]
g[6,y1:y2+1,x2]=( Tw*(w[8]+w[6]) )-g[8,y1:y2+1,x2]
# top Boundary -
g[4,y2,x1:x2+1]=( Tw*(w[4]+w[2]) )-g[2,y2,x1:x2+1]
g[7,y2,x1:x2+1]=( Tw*(w[7]+w[5]) )-g[5,y2,x1:x2+1]
g[8,y2,x1:x2+1]=( Tw*(w[8]+w[6]) )-g[6,y2,x1:x2+1]
# bottom Boundary -
g[2,y1,x1:x2+1]=( Tw*(w[4]+w[2]) )-g[4,y1,x1:x2+1]
g[5,y1,x1:x2+1]=( Tw*(w[7]+w[5]) )-g[7,y1,x1:x2+1]
g[6,y1,x1:x2+1]=( Tw*(w[8]+w[6]) )-g[8,y1,x1:x2+1]
#
##================================
#
for i in range(1,nx):
for j in range(1,ny):
sumt=0.0
for l in range (k):
sumt=sumt+g[l,j,i]
T[j,i]=sumt
T[0:ny+1,0]=0.0
T[0:ny+1,nx]=0.
T[0,1:nx]=T[1,1:nx]
T[ny,1:nx]=T[ny-1,1:nx]
return
In [110]:
def plot (Ra):
strf[0,0]=0.0
for i in range(nx+1):
rhoav=0.5*(rho[0,i-1]+rho[0,i])
if(i>0) :
strf[0,i]=strf[0,i-1]-rhoav*0.5*(v[0,i-1]+v[0,i] )
for j in range(1,ny+1):
rhom=0.5*(rho[j,i]+rho[j-1,i])
strf[j,i]=strf[j-1,i]+rhom*0.5*(u[j,i]+u[j-1,i])
x=numpy.linspace(0,nx,nx+1)
y=numpy.linspace(0,ny,ny+1)
print '************* '
print 'Data for Ra=',Ra
print '************* '
print ' '
plt.figure(figsize=(30,5))
ax = plt.subplot(151)
plt.contourf(x,y,T,30)
title('Temperature Contour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(152)
plt.contour(x,y,strf,30)
title('Streamline',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(153)
plt.contour(x,y,T,30)
title('Isotherm',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(154)
plt.contourf(x,y,u,30)
title('U contour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(155)
plt.contourf(x,y,v,30)
title('Vcontour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()
return
In [111]:
function (1e3,ny/4)
plot(1e3)
function (1e3,ny/2)
plot(1e3)
function (1e3,3*ny/4)
plot(1e3)
In [112]:
function (1e4,ny/4)
plot(1e4)
function (1e4,ny/2)
plot(1e4)
function (1e4,3*ny/4)
plot(1e4)
In [113]:
function (1e5,ny/4)
plot(1e5)
function (1e5,ny/2)
plot(1e5)
function (1e5,3*ny/4)
plot(1e5)
In [113]: